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ABSTRACT 

We present an algorithm that computes the multipole coefficients of the galaxy three- 
point correlation function (3PCF) without explicitly considering triplets of galaxies. 
Rather, centering on each galaxy in the survey, it expands the radially-binned density 
field in spherical harmonics and combines these to form the multipoles without ever 
requiring the relative angle between a pair about the central. This approach scales with 
number and number density in the same way as the two-point correlation function, 
allowing runtimes that are comparable, and 500 times faster than a naive triplet count. 
It is exact in angle and easily handles edge correction. We demonstrate the algorithm 
on the LasDamas SDSS-DR7 mock catalogs, computing an edge corrected 3PCF out to 
90 Mpc//i in under an hour on modest computing resources. We expect this algorithm 
will render it possible to obtain the large-scale 3PCF for upcoming surveys such as 
Euclid, LSST, and DESI. 

Key words: cosmology: large-scale structure of Universe, methods: data analysis, 
statistical 


1 INTRODUCTION 

In the current picture of structure formation, inflation ends 
in reheating, which produces Gaussian random field den¬ 
sity fluctuations in the radiation, matter, and dark matter. 
As a Gaussian random field, the density is described com¬ 
pletely by its mean and 2-point correlation function (2PGF), 
which measures the probability of finding a certain value 
of the density at one point given the density at another. 
However, the subsequent evolution of the density field in¬ 
troduces additional correlations as gravity drives the con¬ 
vergence of overdense regions towards each other. In partic¬ 
ular, a 3-point correlation function (3PGF) is produced by 
this evolution (Bernardeau et al. 2002 or Szapudi 2005 for 
reviews). Since the evolution is itself sensitive to the cosmo¬ 
logical parameters, measuring the 3PGF of galaxies offers an 
independent probe of these parameters. It is typically used 
to break the degeneracy between galaxy bias (encoding the 
fact that galaxies do not trace the matter density field with 
perfect fidelity) and the clustering on a given scale (e.g. as) 
(Gaztahaga & Frieman 1994; Jing & Borner 2004; Guo et 
al. 2014). The 3PGF measurements also can probe primor¬ 
dial non-gaussianity (Desjacques & Seljak 2010); while the 
constraints on this are currently dominated by GMB exper- 
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iments such as Planck, it is expected that the increasing 
quality and number of galaxy redshift surveys will provide 
interesting independent information. 

Since the first measurement by Peebles & Groth (1977), 
numerous studies have presented 3PGF measurements, sum¬ 
marized in Kayo et al. (2004), McBride et al. (2011a, b), 
Guo et al. (2014) and references therein. In this work, we 
present a new algorithm for measuring the 3PGF of galaxies 
through its multipole moments. This decomposition of the 
3PGF was first advanced in Szapudi (2004) and to a limited 
extent (measurement of the monopole moment) used in Pan 
& Szapudi (2005) on Two-degree-Field Galaxy Redshift Sur¬ 
vey (2dFGRS) data. Slepian & Eisenstein (2015) (hereafter 
SE15) found this decomposition to be particularly useful for 
distinguishing linear and non-linear bias as well as isolating 
a possible relative velocity bias. 

Gurrent algorithms, such as that used for the McBride 
et al. (2011) measurement (presented in Moore et al. 2001, 
Gray et al. 2004, Nichol et al. 2006, and Gardner et al. 2007) 
fundamentally scale as the number of possible triangles in 
a survey. If one wishes to measure the 3PGF out to some 
scale Rmax, there are relevant triangles, where 

N is the number of objects in the survey, n is the survey 
number density and = (4/3)7rRmax- The algorithm 

presented in the series of references above, whose most re¬ 
cent incarnation is developed in March (2013), uses multiple 
mrkd-trees. Here “mr” means the kd-tree caches additional 
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information, in this case the number of galaxies within each 
node of the tree as well as the bounding box of the node. 
This algorithm is faster than simply counting all triangles. 
It is particularly effective if the galaxies are close to each 
other, so that there are many triangles whose side lengths 
fall within a given combination of radial bins. 

However, typical galaxy surveys are sparse, particularly 
those mapping the largest volumes. For example, the Baryon 
Oscillation Spectroscopic Survey (BOSS) has an average 
separation of 13 Mpc//i, too large to permit many galax¬ 
ies to be in the same bin. This means the algorithm will not 
be as fast for such large-scale measurements. The use case 
tested in March (2013) is triangles with three sides of 8 Mpc 
each, much smaller than the scales that are well-described 
by linear perturbation theory and hence most useful for cos¬ 
mology. Furthermore, even with the speed-ups coming from 
the multi-tree structure of the algorithm, it is still funda¬ 
mentally scaling as the number of galaxies in the survey 
times the square of the number within i^max (March 2013, 
Figure 21). 

In this paper, we present an algorithm that does better: 
it scales as the number of galaxies in the survey times the 
number within i^max, and so by construction is significantly 
faster than any previous algorithm that is exact in angle. In 
brief, we write the opening angle dependence of the trian¬ 
gles about a given vertex in terms of Legendre polynomials 
of fi -7^2, where these are two unit vectors describing two 
triangle sides. The dot product seems to require explicitly 
considering all pairs of galaxies about a given vertex (i.e. 
third galaxy), but using the spherical harmonic addition 
theorem, this representation can be factored into a prod¬ 
uct of spherical harmonics each depending on only one unit 
vector. Therefore from the spherical harmonic expansion of 
the radially binned density field one can obtain the multi¬ 
pole moments without ever needing to consider pairs about 
a given vertex. This is the central insight of this paper. 

In Section!^ we present the algorithm in more detail, 
and show in Section how this framework goes through 
to the projected 3PCF. Section discusses edge correc¬ 
tion, while Section describes our implementation. Section 
[^computes the covariance of this multipole decomposition 
in the Gaussian random field limit, and Section presents 
the results of using the algorithm on the LasDamas SDSS-II 
Data Release 7 (SDSS-DR7) Luminous Red Galaxy mock 
catalogs. We conclude in Section 


2 THE ALGORITHM 
2.1 Legendre basis 

In this paper, we parametrize triangle configurations by two 
side lengths, ri and r 2 , and the angle between them with 
cosine fi -72. We will decompose the 3PGF as a function of 
these three variables into a sum over Legendre polynomials 
for the angular dependence times radial coefficients encoding 
the side length dependence, as 

C(^i,^2;ri • r2) = ^0(^1,^2)Pz(ri • r2). (1) 

Szapudi (2004) first advanced this decomposition, and he 
puts a factor of {21 + l)/(47r) in front of his analogous ex¬ 
pansion coefficients; we absorb this into Q. 


There are three major advantages to this decomposi¬ 
tion. First, the shape of the 3PGF for fixed side lengths as 
a function of angle is smooth and slowly varying (see e.g. 
Bernardeau 2002, Figure 11), without much fine structure. 
Thus we expect that only a few multipoles will be required 
to capture the angle dependence. Second, this decomposi¬ 
tion provides a natural way to visualize the 3PGF for all 
triangle configurations; one can make several panels for dif¬ 
ferent /, each with all ri and r 2 and amplitudes indicated by 
a color bar, as in SE15. In contrast to many previous works, 
this allows immediate appraisal of the information in all tri¬ 
angles and not just a particular set of configurations (e.g. 
isosceles, two-to-one, etc.) 

Third, as we will see, the multipole moments of the 
3PGF can be obtained with much greater speed than other 
decompositions of the 3PGF. However, in contrast to other 
fast methods, such as tree methods that fix a critical an¬ 
gular scale below which they are approximate (e.g. Zhang 
& Pen 2005) or Fourier methods that choose a grid with 
some minimum spacing, we do not sacrifice accuracy to ob¬ 
tain this speed. Our method is exact in angle. We will bin 
in side length, but even were speed of no concern this would 
be necessary to keep the covariance matrix to a reasonable 
size. 


2.2 Rotation and translation averaging 

The 3PGF describes the number of triangles of a given con¬ 
figuration whose vertices are the galaxies in a survey. While 
nine coordinates are required to completely describe any in¬ 
dividual triangle connecting three galaxies, the 3PGF av¬ 
erages over both translations and rotations of the triangle 
configuration. The presumed losslessness of this averaging 
corresponds to the two usual cosmological assumptions of 
isotropy (rot at ion-invariance about a given point) and ho¬ 
mogeneity (translation-invariance). This ultimately reduces 
the 3PGF to a function of only three variables; as indicated 
already, we will use two triangle sides and the angle between 
them. We will now show explicitly how to go from nine co¬ 
ordinates to three. 

We begin with averaging over rotations. We will show 
explicitly that Legendre polynomials are an angular basis 
for the 3PGF after this averaging. To do so, we first step 
back and write an estimate (denoted by a hat) of the 3PGF 
for a triangle with sides fq, fq extending from a vertex whose 
absolute position within the survey is s. We have 

Im I'm' 

We now wish to average over all rotations of the triangle 
about s. Writing a rotation as R (simply a matrix involving 
the three Euler angles), we have 

Ciso(ri,r2;fi •r2;s) = (ri,r2;s) 

Im I'm' 

X J dRYim{Rf\)Yi1^,{Rf2), (3) 

where subscript “iso” abbreviates “isotropy.” Noting that 
Yim{Rr) = ^mM^ZM(r), where D\^m is a Wigner ma- 
trix (e.g. Arfken, Weber & Harris 2013 (hereafter AWH13), 
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equation (16.52)), we find 

6so(ri,r2;ri •f2;s) = {ri,r2-,s) 

Im I'm' 

X (^ 2 ) f dR M'■ (4) 

MM' ^ 

The integral over Wigner matrices is simply evaluated by 
orthogonality (e.g. Brink & Satchler 1993, Appendix V) as 
87 r^/( 2 / + , 6^ the Kronecker delta. Using 

the spherical harmonic addition theorem (AWH13, equation 
(16.57)) 

a(ri • r2) = y] Yi^ih)YUt2), ( 5 ) 

m = — l 

and defining 

Ci{ri,r2-,ri •■r2;s) = 27r C/T’^Cn, r2; s) (6) 

m 

we find 

6 so(ri,r 2 ;ri • r 2 ]s) = '^Ci{ri,r 2 ] s)Pi(ri • f 2 ). (7) 

i 

In what follows we drop the subscript “iso” as we will always 
be considering the isotropic 3PCF. 

We now move to averaging over translations. Recalling 
that s is the vertex of the triangle from which the two sides 
given by 75,72 extend, the densities on a particular triangle 
of points will be S{^S{fi + s)S{f 2 + ■^. Averaging over trans¬ 
lations means allowing every point in the survey to serve as 
the vertex s, so we must integrate over d^s. We thus find 
that the radial coefficient of the 3PCF is 

Cl{ri,r2) = ^ J (fsCi{ri,r2;s), (8) 

where V is the survey volume. 


2.3 Radial binning 

Our algorithm will bin radially (denoted with a bar), so we 
seek 

Ci{ri,r2) = j r‘^r‘^drdrCi{r,r)^(r]ri)^{r]r2), (9) 

with $ a binning function demanding that we are in the bin 
given by its second argument. Binning averages the radial 
coefficient over some interval in each side length, and in 
that sense is not lossless. It is also necessary for the speed 
advantage of our algorithm, as will become clear shortly. 

We will not compute using equation Rather, we will 
bin radially around each possible origin s before averaging 
over rotations and translations, so it will be useful also to 
define the binned estimator before translation-averaging as 

C;(T’i,r 2 ;s) = 

J dQidQ2S{s)5{ri-,fi-,^S{r2-,f2-,^Pi{fi ■ f2), ( 10 ) 

where 

= J r‘^dr ^(r;ri)6{fi s) (11) 

is the radially binned density field about an origin s. 
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Hence in practice we never compute (i (r, r ') u sing equa¬ 
tion 1 ^, but rather measure Q via equation ( [l^ and then 
compute 

Ci{ri,r 2 ) = y Jd^s Ci{ri,r 2 -,s) ( 12 ) 

as the radially binned multipole coefficients of the 3PCF. 


2.4 Accelerating with spherical harmonics 


A direct way to measure Q would be to sit on every possi¬ 
ble origin and compute the angle between pairs of vectors 
pointing to all possible sets of two galaxies out to the ra¬ 
dius Rmax to which one wishes to measure the 3PCF. This 
scales as N{nVR^^^y. As discussed in the Introduction, this 
scaling applies to other algorithms as well (e.g. the Gardner 
(2007) and March (2013) kd-tree approach), fundamentally 
because the number of possible triangles within Rmax with 
one vertex fixed scales as (nRmax)^. 

However, as is the case for angular power spectra, we 
can exploit a property of multipole decompositions to enor¬ 
mously accelerate the measurement. We can use the spheri¬ 
cal harmonic addition theorem 0 to decompose the Legen¬ 
dre polynomial into factors that depend only on one angular 


variable each. Inserting this into equation (10), we find 


6(n,r2;s) = ^(5(s) ^ J dQi S{ri-,fi-,^Yimiri) 

m= — l 


j dO.2 S{r2-,r2-,s)YMf2). 


(13) 


This equation immediately shows how to reduce the 
quadratic scaling in the number density to a linear scaling. 
The two angular integrals have now been separated, and 
each simply asks for a particular expansion coefficient of the 
density field (as a function of angle alone) in spherical har¬ 
monics, in a fixed radial bin. In other words, if we compute 
for each radial bin r 

aim{r;s)= JdQ 5{r;f-,s)YMr) 

= J dQYi^{r) J r'^drA{r';r)6{r + s) (14) 

we can construct all combinations dictated by ri and r 2 , 
without ever needing to do an (!l(n^) operation. Explicitly, 
inserting equation (14) into equation (13), we find 


^ 1 ^ 

Ci(»'i>'r2;s) = ^<5(s) y] ai^(ri;s)a*^(r2;s). (15) 


This is why radial binning is essential for the speed-up of 
our algorithm; we can precompute the aim{r;^s in each 
radial bin. For Abins, we then only need to construct (AbinsT 
l)A’bins/2 combinations of these coefficients. A schematic 
about a single possible origin is shown in Figure 

For a 3PCF measurement, one might use a bin width ~ 
10 Mpc, and so if one measures out to 200 Mpc there will be 
only 210 distinct bin combinations. Meanwhile, computing 
the aim{r;^s themselves takes only as long as performing 
the integral which should scale as . 

We still must integrate over all possible choices of ori¬ 
gin as dictated by equation Because the galaxies are 
discrete, this will reduce to a sum with N terms. Thus our 
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Figure 1. Our algorithm sits on each galaxy in the survey, here 
marked with a white X, and computes the spherical harmonic ex¬ 
pansion of the density field in concentric spherical shells (radial 
bins) around that point via equation in- The ai^s can be com¬ 
bined to yield the multipole moments around this galaxy (sum 
over m, equation ( |15| )) and then translation-averaged to yield Q 
for the survey. 


algorithm will scale as N{nVR^^^): linear in both the total 
number of galaxies and the number within a sphere of radius 
i^max, and a factor of order faster than the naive 

counting approach. Our algorithm thus provides a route to 
the 3PCF that on large scales is no more computationally in¬ 
tensive than calculating the multipole moments (standardly 
calculated are monopole and quadrupole) of the 2-point cor¬ 
relation function (2PCF). 

Finally, to obtain the spherical harmonic coefficients 
of the galaxy density as in equation one might think 

a spherical harmonic transform is required. This scales as 
, Ng the number of spatial grid cells on the surface of 
a sphere. A large number of grid cells is necessary for ac¬ 
curacy even if there are very few galaxies, much as a small 
Ak is needed when taking a numerical Fourier transform to 
avoid ringing. However, because only low-order multipoles 
are needed here (I < 10), we can avoid this transform and in¬ 
stead directly evaluate the a^mS, which are simply spherical 
harmonics evaluated at angles given by a galaxy’s location 
with respect to a given choice of origin. The required YimS 
can be easily computed using the Cartesian expressions for 
the spherical harmonics (e.g. AWH13, equation (15.139) and 
Table 15.4). Indeed, about a given origin, the Cartesian com¬ 
ponents x/r, y/r, and zjr and their powers for each galaxy 
can be pre-calculated just once and subsequently combined 
to form all of the required multipoles. 


3 PROJECTED 3PCF 

Redshift-space distortions (RSD) are differences between the 
true position of a galaxy along the line of sight and its posi¬ 
tion as inferred from assuming its redshift is purely cosmo¬ 
logical. They arise from peculiar velocities, ultimately gener¬ 
ated by the growth of large-scale structure (Hamilton 1998, 
for a review). The projected 3PCF is insensitive to these 
distortions because it is integrated along the line of sight. 
Below we show how our approach extends to measuring it. 

We work in the flat-sky approximation, where there is 
a single line of sight to all galaxies in the survey. Sitting 


around a given central galaxy and projecting corresponds 
to drawing cylindrical shells around that central with bases 
that are concentric annuli. All of the galaxies in a given 
cylinder project down into the cylinder’s base annulus. We 
thus have a planar problem with circular symmetry. 

This permits simplification of our spherical harmonic 
basis. Recall that 

where here cj) are the angular coordinates of a galaxy in 
the system where the central is at the origin. Since all the 
(projected) positions are coplanar with the central, the sep¬ 
aration along the z-axis is zero, so cosO = 0. Defining 


^Im — 


2/ + 1 (/ — m)! 

Att {Im)\ 




(17) 


we see from equation (13) that the muitipoie moments of the 
projected, radiaiiy binned 3PCF wili simpiy invoive Fourier 
coefficients of the projected, radiaiiy binned density fieid 
weighted by bim'- 


0 ,proj (^ 1 , ^2; ^ 


E 


#i(^(ri; 0i; s)F 


im4>i 


J 


d(t)25{r2] 02 ; 


— im4>2 


(18) 


Above, the integrais over 0i and O 2 of equation (13) have ai- 
ready been performed using that the projected density field 
is only non-zero at 0 = 7r/2. 

With this in mind, we observe that if one is solely inter¬ 
ested in the projected 3PCF, it is probably optimal simply 
to use the Fourier basis directly. One parametrizes the pro¬ 
jected 3PCF estimator about a given central as 

4roj(ri,r2;6>i2;s) = X ^ 2 ; (19) 


and writes the exponential as 

^ (e*"*®!)* ( 20 ) 

where 0i and O 2 are now angles in the plane in polar co¬ 
ordinates, with O 12 = O 2 — Oi. Using orthogonality of the 
plane waves, one may then extract the expansion coefficients 
Cproj,m(ri,r 2 ; ^ in equation (19) as 


Cproj,m (^15 *^ 2 ; — 


(27r) 


■/ 


dOi S{ri; Oi; s)e^' 


/ 


d6i2 <5(r2;6i2;s)e-*’"®E ( 21 ) 


Just as in the non-projected case, these integrals can be 
explicitly evaluated using the Cartesian expressions for the 
exponentials, and precomputing x/r and y/r. Again, one 
never explicitly considers pairs of galaxies about a given 
central; one simply constructs the coefficients 


Cm(r’;s) = J dO 

for all radial bins, then computes 

4roj,Tn(ri,r2;s) = 5{s)Cm{ri\s)c,n{r2s)l{2'K)'^ (23) 

for all desired bin combinations, and finally averages over 


( 22 ) 
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translations by integrating out s. We should note that Chen 
& Szapudi (2005) advanced a similar scheme to measure the 
3PCF of Cosmic Microwave Background (CMB) maps, anal¬ 
ogous to the projected 3PCF since both are on 2-D mani¬ 
folds. However their method evaluates the Fourier transform 
of the (continuous) temperature anisotropy map by grid- 
ding, whereas here we suggest the (discrete) galaxy density 
field be Fourier-transformed using direct evaluation of the 
Cartesian expressions for x/r and y/r. 


4 EDGE CORRECTION 


Surveys have jagged and complicated boundaries, and these 
can produce a spurious contribution to the 3PCF that is 
the signature of the survey geometry rather than physics 
the survey hopes to probe. This spurious contribution must 
be removed. In Fourier space, boundaries lead to Gibbs phe¬ 
nomenon ringing in the bispectrum, and are challenging to 
remove. However, in configuration space, edge correction is 
fairly straightforward for popular estimators (see Kayo et al. 
2004, Appendix, for comparison of several). 

We focus here on the Szapudi & Szalay (1998) estima¬ 
tor, which Kayo et al. (2004) find preferable to the others 
they consider; it has now become the standard in the field. 
It is 


NNN 

RRR 


(24) 


with N = (D — R), D the data and R the random counts. 
Note that, if one inserted N/R for S in Section]^ one would 
need to compute integrals of this fraction against the spher¬ 
ical harmonics, requiring definition of N/R at every point 
in space. However, the estimator ( [24| really represents the 
function 


Craw(Fi,r2,r3) 


N{fi)N{r2)N(r3) 

R{fi)R{f2)R{f3) 


(25) 


averaged over rotations and translations with weights w = 
R{fi)R{f 2 )R{r 3 ) 0 , which in the shot noise limit is just in¬ 
verse variance weighting (we include radial binning repre¬ 
sented by 0). In short. 




/ d^fi(ff2d^f3 re(ri,r2,r3)Craw(TL,r2,r3) 

f d^fid^f2d^f3 w{fi,f2,r3) 

J d^rid^r 2 d^ f 3 ONNN 
J d^fId^f2d^r3 ORRR 


(26) 


Thus the estimator (24) should be interpreted as demanding 


the triple count NNN divided by the triple count RRR. 
Therefore we can insert N and R separately in turn for S 
in Section!^ processing random and data catalogs serially. 
The division required can be done as a post-processing step. 
We now turn to how this division translates to the Legendre 
basis. 


4.1 Edge correction in the Legendre basis 

Working now in our Legendre basis, we have 

C = '^Ci{ri,r2)Piiri-f2) (27) 

I 
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NNN = E J^j{ri,r2)Pj(fi ■f2) 


( 28 ) 


and 


RRR = ^ Til/ (n, r2)Pv (r. 


1 -^ 2 ) 


(29) 


Inserting the multipole expansions (27)-(29) into the esti¬ 


mator ([24| and multiplying through by RRR we find 

Y,Pi<iPviA ■ r2)Pi{ri ■ f2) = VP,(ri • ra). (30) 


Using a linearization formula for the product of two 
Legendre polynomials (Ferrers (1877), Adams (1878), Neu¬ 
mann (1878), Park & Kim (2006); SE15 equation (All)) we 
find, with angular arguments suppressed. 


^77K,(2/ + l) 

I'lj' 


0 


J 

0 


Pj, = J2-^jPr ( 31 ) 


The Wigner 3j-symbol above describes angular momentum 
coupling; see e.g. Brink & Satchler (1993) or AWH13. The 
vector addition of angular momenta means that the upper 
row must satisfy triangle inequalities, so |/ — /^| ^ j' ^ I 1' 
and at fixed I and I' the sum is hnite. Using orthogonality, 
separating out the = 0 term, dividing through by 77.o, and 
defining fif = IZif /IZq, we obtain 

f = a + Eo(2‘+i)E(; Sp"' 

I i'>o ^ ^ 


For a boundary-free survey the random field would gen¬ 
erate only a monopole (77.o), leaving only (/k on the right- 
hand side; this is the limit where there is no need for edge- 
correction, but just division by the randoms. The form of 
equation (32) suggests that this problem can be cast as a 


matrix multiplication, so we define the multipole coupling 
matrix AT with elements 

M,t = {2k+l)J2{ 0 0 o) 

l'>0 ^ ' 


Note that while these matrix elements describe the off- 
diagonal couplings of different multipoles to each other, they 
need not be zero along the diagonal. A given multipole in 
the data may couple to that same multipole in ^ because 
the 3j-symbol allows I = k for /' > 0. But the dominant cou¬ 
pling of a given multipole in the data to the same multipole 
in is described by in equation (32), since the fi/ are 


expected to be much less than unity. This term translates 


to the identity matrix I. The edge-correction equation (32) 
thus becomes 


A7■/7^o = (/ + M)C = AC 


(34) 


where A/ = (A/o, A/i, • • • , A/z^^x) analogously for C The 

system of equations this represents can then be solved for ( 
by matrix inversion. 

To explore this matrix for a realistic use case, we use 
the LasDamas SDSS DR7 real space mock catalogs, using 
15 radial bins and a maximum scale of 90 Mpc/h (further 
details are given in Section 1^ We show Mki for a partic¬ 
ular bin in (ri,r 2 ) in Figure]^ and the leading order edge 
correction factor fi in Figure^ AT is not symmetric, but 
Mki/{2k + 1) is; this is why the upper off-diagonal, where 
k > I, exceeds the lower in Figure 
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4.2 Solving the edge correction equation 

There are two approximations implicit in our approach to 


solving equation (34). First, to obtain a given matrix element 


Mki, formally one requires fi^ for all values of I'. However, 
for I' > 1 these factors fall rapidly. For the LasDamas real 
space mock catalogs for which we present results here, they 
are < 0.1% by I' = 10 even for the largest-scale radial bin 
combination (the values are listed in the caption to Figure 
1^, so we simply truncate the series there. If one wished one 
could easily expand our code to measure higher multipoles of 
the randoms at the cost of slightly more computation time. 
However we expect that going to = 10 will already render 
the edge correction error negligible compared to the total 
error budget. 

Importantly, the smallness of the fi/ for I' > 1 means 
that the coupling between multipoles k and I is nearly diag¬ 
onal. Coupling between multipoles separated by more than 
one angular momentum step is suppressed as /2 or higher 


because the 3j-symbol in the coupling matrix elements (33) 
requires that I' > \l — k\. 

The second approximation relates to the matrix inver¬ 


sion when we solve equation (34). Formally one has an in¬ 


finite dimensional matrix where at fixed /, all k enter the 
correction. Thus this matrix will not be square (and hence 
invertible) unless we go to an infinite number of I as well. 
However, in practice the matrix is so diagonally-dominant 
that we believe it is accurate enough simply to invert the 
sub-matrix given by truncating I and k at some maximum 
multipole. We verify this approximation by constructing Mki 
using solely the dominant fi edge-correction factor, letting 
k and I go to 2/max, inverting, and comparing to the result 
where both go to /max- Were the matrix purely diagonal, 
truncation would not affect the inverse at all. In the limit 
where only fi is non-zero (in reality, it does dominate the 
other edge correction factors), the matrix is tridiagonal, and 
so truncation at /max affects Omax order /i, Omax-i 
der /f, and at order 


4.3 A model for the edge correction factors 

Using a simple toy model, we can estimate the edge correc¬ 
tion factors fi to confirm that they really should be small. 
Consider a spherical ball of random galaxies with radius R 
about a given central, and assume this sphere is cut by a 
planar survey boundary. Orient the z-axis perpendicular to 
this boundary, with the central galaxy a distance z from it. 
The problem now has symmetry about this axis, so we need 
only compute the m = 0 spherical harmonic coefficients aim] 
Yio = with fi = cos^. For a galaxy at 

distance R from the central, there will be some critical angle 
with cosine fic = z/R such that, for smaller /x, the galaxy is 
outside the survey. We have 


aio = 27V 


27V 




2 / + 1 
4:7V 


[-R + i(Mc) -R-i(Mc)] • 


(35) 


2 / + 1 

We used the recursion formula (2n + l)Pn{/a) = 
d/dfi [Pn-\-i{/a) — Pn-i{/a)] to evaluate the integral and noted 
that the terms at the lower bound cancel off because they 
have the same parity. We now compute the Q = a^o/(47r) 


required by equations (14) and (13) and average over fic 
(denoted by angle brackets). We have 


(‘^*) = 4(2ZTT)/ 

- 2Pi+i(/ic)Pi-i(Atc) + Pi^_i(/ic)]- (36) 


Since each term above has even parity, we can integrate from 
— 1 to 1, divide by 2, and then invoke orthogonality, to find 
that 


( 




1 

2(2/ + 3)(2/ - !)■ 


(37) 


Finally, we compute 



= 7/12 explicitly, to find that 



/i = 17.14%, /2 = 4.08%, /3 = 1.90%, /4 = 1.11%, faffing 
to /lo = 0.196%. It should be kept in mind that in a 
large survey volume such as SDSS, many centrals will have 
spheres around them that do not impinge on a large-scale 
survey boundary at all, further reducing these factors; for 
instance, for the SDSS BOSS DRIO footprint only of order 
20% of spheres impinge on a boundary, so our rough esti¬ 
mates should be scaled down by a factor of 5. On the other 
hand, the true survey mask is far more complicated than the 
simple planar boundary model above, so this model should 
not be taken too literally. 


5 IMPLEMENTATION 

We next describe our C++ implementation of the ideas in 
Sections and The basic program flow is to loop over 
each central galaxy. For each, we find all neighbors within 
Rmax and accumulate the aim for each of the radial bins. 
Once finished with the neighbor finding, we compute all of 
the bin cross-powers and add them to our accumulators as 
a function of bins ri and r 2 and multipole /. All of the ac¬ 
cumulations include a user-supplied weight per galaxy. 

We accelerate the finding of neighbors by sorting the 
particles into a grid, so that the search for neighbors need 
only consider grid cells that include some point closer than 
Rmax. Ideally one wants the grid spacing to be a few times 
smaller than Rmax, so that the inefficiency of doing a cubic 
search for a spherical region is mild. One also wants the 
grid spacing to be large enough to contain at least several 
particles, so that the overhead of storing and accessing the 
grid is modest. These criteria are not hard to satisfy: for the 
LasDamas mocks, we use a grid spacing of 50 Mpc/h when 
searching out to Rmax = 90 Mpc/h; this typically contains 
a dozen galaxies (and somewhat more random points). For 
the SDSS-HI Baryon Oscillation Spectroscopic Survey, the 
density is three times higher. 

Once a neighbor is found, we need to add its contri¬ 
bution to the spherical harmonics. We do not use an an¬ 
gular binning to compute the spherical harmonics. Rather, 
as mentioned in Section we use the fact that the spher¬ 
ical harmonics can be written as powers of the Cartesian 
coordinates of unit vectors. In particular, for a unit vector 
r = {x,y,z), we can write Yim{r) as a polynomial of terms 
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Figure 2. The leading-order edge correction factor in equation 
( |32| ) for the LasDamas SDSS DR7 real space mock catalogs with 
15 radial bins out to 90 Mpc/h. Higher /;/ = IZi//IZ q fall off very 
rapidly. Even the leading order coefficient is small. This means 
that one does not need to measure many multipoles of the ran¬ 
doms to obtain a highly accurate edge correction: since the higher 
fl/ fall off so rapidly they contribute very little to the matrix M 
that must be inverted (equation (|34|)). As we expect, fi becomes 
larger at larger scales, as larger scale triangles are more likely to 
impinge on a survey boundary. 

of the form where ^ 1. To compute all aim up 

to multipole order p, we therefore accumulate sums over all 
neighbors of the Cartesian powers with ^ p, 

using the unit vectors of the separation of the neighbor from 
the central galaxy. There are {p + l)(p + 2)(p + 3)/6 such 
power combinations for each radial bin. Having finished with 
all neighbors, we convert these powers into the aim using the 
appropriate coefficients from the spherical harmonics, then 
form all of the bin-to-bin cross powers. 

For values of p of order 10, the computation of the 
Cartesian powers is much faster than doing the spherical 
harmonic transform of a fine angular grid. This is particu¬ 
larly true because we use custom assembly code, supplied 
by Marc Metchnik as part of the Abacus project (Metch- 
nik & Pinto, in prep.), to accumulate these powers using 
Advanced Vector Extension (AVX) instructions. In double 
precision, 8 neighbors are computed at once, using two sets 
of AVX registers. 


Figure 3. The multipole coupling matrix elements ( |33| ) at each 
k and I for the largest combination of radial bins we test here. 
This illustrates that all of the couplings are < 10%, even for the 
largest scales we test, which should have the largest correction 
factors as they are most likely to impinge on a survey boundary 
(see Figure]^. While the diagonal appears zero in this plot, it 
is actually just small, as we discuss in the main text. The fi 
entering the matrix elements Mj^i for this radial bin combination 
are fi = 7.18%, /2 = 1.59%, fs = 0.783%, U = 0.428%, /s = 
0.318%, fe = 0.22%, fr = 0.166%, fg = 0.078%, fg = 0.078%, 
and /lo = 0.051%. 
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Though we do not present the 3PCF measurement here, 
we also run our algorithm on the SDSS-III BOSS DRIO data. 
In the North Galactic Cap footprint for the CMASS sam¬ 
ple, we consider the RRR count of 642,619 random particles. 
We count 6.7 billion pairs with Rmax = 200 Mpc//i, an aver¬ 
age of 10,400 neighbors per central, divided into 10 linearly 
spaced radial bins. Using p = 10, the code runs in 170 sec¬ 
onds on a 6-core 4.2 GHz i7-3930K. If we use p = 0, thereby 
reducing the problem to the pair finding and a simple accu¬ 
mulation per radial bin, then the code runs in 53 seconds. 
Loading the particles and sorting them into the grid is a 
small fraction of that total, so we infer that each pair found 
and processed at p = 0 takes about 200 clock cycles. Given 
p = 10, we have 286 powers to track per neighbor, each 
requiring a separate multiply and add. Hence, we are com¬ 
puting about 3.8 trillion double-precision operations in 120 
extra seconds, a rate of 32 double-precision GFLOPS. This 
is about 30% of the maximum performance of the GPU (as¬ 
suming 4 double precision operations in AVX per clock cycle 
per core), a high mark for a practical calculation. The code 
is sustaining 22 GFLOPS for the full problem, including the 
pair finding. 

A two-point correlation function code would only need 
to count half as many pairs, since the particles are indistin¬ 
guishable in that application, and so at these speeds would 
take of order 53/2 = 27 seconds. We therefore find that our 
computation of the three-point correlation function up to 
I — 10 is only about 170/27 6 times slower than the equiv¬ 

alent two-point correlation function calculation. We would 
expect an explicit counting of triples to be about 3,000 times 
slower than the two-point pair counting, given the 10,000 
neighbors (divide by a factor of 3 for the number of indis¬ 
tinguishable triples compared to indistinguishable pairs). As 
our method is only six times slower than a two-point mea¬ 
surement, it is a factor of five hundred faster than an explicit 
triple count for this large-scale example. 

In any method that compares the data points to a ran¬ 
dom set, we have to consider the effect of Poisson noise in the 
randoms. For example, in the Landy-Szalay (1993) estima¬ 
tor for the two-point function, ^ = {DD — 2DR — RR)/RR, 
we will have noise in the data-random (DR) and random- 
random (RR) counts that would go to zero in the limit of 
infinite numbers of random points. One therefore usually 
wants to use many more randoms than data (but note the 
important optimization presented in Padmanabhan et al. 
(2007) in which one fits these counts to smooth functions 
of scale so as to reduce the Poisson noise). A common in¬ 
efficiency, however, is to use the same number of randoms 
for each of the terms. This results in spending far too much 
computational resource on RR, whose Poisson noise will be 
dwarfed by the DR noise. For example, if the number of 
randoms is m times the number of data points, then (as¬ 
suming uniform galaxy weights) the variance on RR will be 
1 /of that of DD since the number of RR pairs is w? the 
number of DD pairs. In contrast, the variance of DR will 
only be reduced to (1/2) (4/m) of that of DD, the factor of 
4 comes from the factor of 2 in the Landy-Szalay estima¬ 
tor, and the 1/2 enters because the DD and RR pairs have 
double the variance since each pair is counted twice. Mean¬ 
while the work in the two terms is scaling as m^/2 and m, 
respectively. 

A simple way to avoid this is to compute the DR and 


RR counts with a smaller set of randoms and then repeat 
this numerous times, averaging over the answers. By choos¬ 
ing the number of randoms in each set, one can optimize 
the work. For example, in the above two-point case, at fixed 
total work, the number of random catalogs A'cat one can 
use scales as (m^/2 + m)~^. The total variance scales as 
the variance per random catalog divided by Neat, so as 
(2/m+ l/m^)(m^/2 + m). This is minimized for m = 1, i.e., 
it is optimal to use random catalogs equal in size to the data 
set. A further advantage of this method is that in addition 
to averaging all of the sets to get the best answer, one can 
compute the variance to explicitly measure the contribution 
of the random catalog density relative to one’s estimate of 
the irreducible on-sky variance. 

For our three-point algorithm, the work scales as 2m + 

while the Poisson variance of (D — R)^ for each random 
catalog scales as 1/m^ for RRR and Im‘^ x (2/6) for DRR 
and DDR; the 3^ enters due to the 3 in the Szapudi-Szalay 
estimator, while the 2/6 comes from a 6-fold counting sym¬ 
metry in DDD and RRR compared to a 2-fold one in DDR 
and DRR. The total variance is thus 3/m + 3/m^ + Ijm^. 
At fixed total work Neat scales as (2m + m^)“^, and so the 
total variance scales as (3/m + 3/m^ + l/m^)(2m + m^). 
This is minimized for m = 1.76 but with only 1% variation 
between m = 1.5 and 2. 

We implement this strategy in our three-point method 
by supplying a single list of particles, with the randoms con¬ 
catenated to the data but with negative weights. Notation- 
ally, this is = D — R, as in Section We then compute 
the three-point correlations of this N list. We then re-run 
repeatedly with new random points R. We avoid the small 
amount of repeated counting of the DD pairs and DDD 
triples by the following trick. We first run the code with 
only the data particle list and save a file that contains the 
Gartesian multipoles for each radial bin and each primary 
particle, in the enumerated order of the particles. When next 
running with D — R lists, whenever a data particle is the pri¬ 
mary (as marked by its having a non-negative weight), we 
initialize the multipole accumulators with the saved values 
and then skip any secondary particles that are also from the 
data list. The resulting sums pass transparently to the rest 
of the analysis code. 

We also run a separate case with only the randoms, so 
that we can compute the denominator and edge-correction 
terms in equation (32). This requires much less precision, as 
the denominator of the estimator is much larger than the 
numerator for large-scale correlations. We therefore do this 
with only a single set of random points. 

Finally, we have also written a Python implementation 
of the algorithm presented here and tested it on a periodic 
box with sides of 400 Mpe/h containing 20,000 galaxies 
(roughly the SDSS BOSS number density). Rather than us¬ 
ing gridding, this code exploits kd-trees for galaxy finding, 
using a fast G implementation (wrapped to python) in the 
spatial library within scipy. We verified the accuracy of this 
code on a sample of 500 galaxies by comparing with a sim¬ 
ple direct-counting algorithm that just counts triplets and 
then projects onto multipoles. This provides an important 
cross check on our spherical harmonics since the simple triple 
counting never uses spherical harmonics. We then ran both 
the multipole Python code and the multipole G++ code on 
a larger, 20,000 galaxy sample to verify the G++ code. Run- 
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time for the Python version on a dual core (2014) MacBook 
Air was about 30 minutes; since the box is periodic, scaling 
to larger numbers of galaxies is linear. 

6 COVARIANCE MATRIX 

Parameter fitting requires weighting the data points accord¬ 
ing to how independent they are, with two highly indepen¬ 
dent points contributing more than two less independent 
points all else equal. The covariance matrix describes how 
independent the measured multipoles at each (ri,r 2 ) are. 
For our algorithm to be useful, we must show that the co- 
variance matrix can be controlled; here we compute it with 
this end in mind. The general 3PCF covariance has been 
computed before (Szapudi 2001) as a 6-D integral, but it 
is not straightforward to obtain the covariance of our mul- 
tipole decomposition from this result. Here we derive the 
covariance for our multipole decomposition and show that 
it can be reduced to a sum of 2-D integrals. This reduction 
offers a significant improvement in the computation speed 
possible at a given accuracy. 


6.2 Full covariance 

We now obtain the covariance of our multipole decomposi¬ 
tion of the 3PCF[^Here we begin with an estimator for the 
translation-averaged but not rot at ion-averaged full 3PCF; 
we will project onto multipoles (which also averages over 
rotations) and bin radially later. 

C(n,^ 2 ) = J ^ (5(s)(5(s + ri)(5(s-I-r 2 ) (45) 

For a Gaussian random field, < ^ >= 0- The covariance is 
thus 

< C(ri,r-' 2 )C(r 1 ,r 1 ) >= J 

X ^ exp [ - z(fc • s + (s + ri)+p • (s + r 2 ) 

kqp,k' q' p' 

k' ■ s q ‘ {s Vi) p ■ {s f^)] 

X (6{k)6{^6{]^6{k')6{q)6{p)'^ . (46) 

Peforming the integrals over and we have 


6.1 Conventions 


< C(ri,f2)C(ri 


kqp,k' q'p' 


We begin with some definitions and conventions. While we 
wish to compute the covariance matrix of the configuration 
space 3PCF, we will end up working in Fourier space to do 
the computation because simplifications are available there 
by appeal to the power spectrum. We define the Fourier 
transform as 

5(fc) = J cTrd{r)A^'^ (39) 

with inverse 

= I Wf \ ( 40 ) 

For the earlier stages of our computation we will in fact need 
to use the discrete Fourier transform and its inverse, defined 
as 

5{k) = ^5{r)e’‘'' (41) 

r 

and 


X exp [ - i{q- fq + p • r2 + • fq + ^ • (^)] 

X (^S{k)mm5{P)S{^)S{^)) (47) 

where is a Kronecker delta whose argument is the sum 
of the subscripted vectors. We now use Wick’s theorem to 
reduce the 6-point expectation value to triple products of 
2-point functions; this is where Gaussianity enters. We need 
to consider all possible contractions. 


C(Ti,r2)C(fq,r2) 


_ cK rK -i[q-ri+p-r2+q -Vi+P -^2] 

ye 2-^ ^kqp^k'q'p'^ 


kqp,k'q' p' 


X |(w')(PP')(fcfc') + {pq){QP){kk') + {kq){qk'){pp) 

+ {kp){qk'){pq) + {kq){pk'){qp) + {kp){pk'){qq)'^ 

(48) 


Sip = y^S{k)t 


ik-r 


where these discrete transforms are over a volume V — 
with quantized wavenumbers such that k 

We define the power spectrum as 

P{k)d^{k + k') = ^(^dik)5{k'))-, 

6^ is the Kronecker Delta, unity when its argument is zero 
and zero otherwise. One can check easily that this definition 
allows one to recover the familiar relation that the correla¬ 
tion function is the Fourier transform of the power spectrum. 
We will also use the fact that 


J (frA^^~^^'k^ = VS^{k + k'). 


(44) 


Finally, note that one can convert from the discrete to 
the continous case by replacing (1/K) '^^k with f d^k/{27r)^. 


(42) where parentheses represent contractions of 6 evaluated at 
the arguments in the parentheses. Using equation (43), the 
term in curly brackets above becomes 


27rn/L, 

{■■} 

(43) 

+ ^kq'^ 


= P{q)P{p)P{k)V^ 


cK rK cK . J-K J-K cK 
^qq' ^pp' ^kk' > ^pq' ^qp' ^kk' 


cK riC cK 
^kp' ^qk' ^pq' 


cK riC cK 
^kq' ^pk' ^qp' 


cK riC rK 
^kp' ^pk' ^qq' 


(49) 


^ The techniques used here can also be used to compute the 
covariance of the full 3PCF, without projection onto the Legendre 
polynomials, in terms of 2-D integrals, but one obtains an infinite 
sum over angular momenta. We have not assessed what error 
truncating this sum might induce. There may be applications 
where this 2-D integral representation, despite the infinite sum 
over angular momenta, could be preferable to the 6-D integral 
expression of Szapudi (2001). 
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Doing the sums over k', q', and p' in equation (48) we find 


C(n,r2)C(rq,r^2)) = ^ P(g)P(p)P(fc) 

kqp \ J 


-i[q-ri+'p-r2\ 


X < e 




(50) 


Notice each pair of exponentials in the curly brackets is obvi¬ 
ously symmetric under switching ry o r^. Also notice from 


the first line that equation (50) is symmetric under fq ^ fq 
if we also flip q and p. Applying this to all of the terms in 
curly brackets too, we find 6 + a + e + / + c + d if we had 
originally labeled each exponential in the curly brackets as 
a + 5 + c + (i + e + /. Hence the equation has the desired 
symmetries. 

Converting this into an integral we have 

Cov = (c{ri,r2)Q{riy2)) = ^ j ^ ^2 ^9 ^ P{p)P{<l)P{k) 
X (2-7r)^ 5^' (q + P + k^ g-*[9 >-i+p »-2] I • • • |, (51) 


where above we have not rewritten the terms in curly brack¬ 
ets from equation (50). 


6.3 Projection onto Legendre polynomials 

We now consider the covariance projected onto multipoles, 
defining 


Covw{ri,r2;r'i,r2) = 


(2/ -|- 1)(2/^ + 1) 


/ 


(47r)4 

X Cov(fi,f2;fq,fq)^z(^i • r2)Pi'{K • ^ 2 )- 


(52) 


Noticing that in equation ( |51| ) the exponentials contain the 
only f dependence, we first define the projection of one ex¬ 
ponential onto one Legendre polynomial as 


i,z(fq,/ci;f2,/c2) 


(47r) J 

= (2/ + i)(-i)*j;(fci,fc2)Pi(fci -fo). 


(53) 


We will have this factor from projecting the exponential out¬ 
side the curly brackets in equation ( |51[ ) , and then six analo¬ 
gous factors within the curly brackets from projecting each 
exponential of onto Pvir'i • r^)- 

We have defined J'i{x,y) = ji{xri)ji{yr 2 ) and will also 
use Ji', {x, y) = jV (xri)jV (yr^). We performed the projection 
integral by expanding the exponential in spherical harmon¬ 
ics using AWH13 equation (16.61) and expanding the Leg¬ 
endre polynomial in spherical harmonics using the spherical 
harmonic addition theorem the integral can then be 
evaluated by orthogonality. 

Writing out the projection integrals explicitly using 


equation (53), we thus have the projected covariance as 


Goyw = ^ j ^ ^ Pip)Pi<l)Pik) {2Trf + q + k) 

X {21 + 1){21' + l){-iy+^' Ji{q,p)Pi{q ■ p) 


X ^Jv{(l,P)Pi' (9 • p) + Jv {p, g)Pi' (p • 4 ) + Ji' {k,p)Pi' {k ■ p) 

+ Jl, {p, k)Pi, {p-k) + J{, {k, q)Pi, {k-q)p Jl, {q, k)Pi, (^ ^ |. 

(54) 


Note that in equation (54) the Legendre polynomial depen¬ 
dence is the same for each of the first pair in the curly 
brackets, the second pair, and the third pair because the 
dot product is symmetric. Thus we have three possible an¬ 
gular integrals to do, corresponding to these three pairs: 


C™ = / d%dn,dQtPliq ■ p)Pv {q ■ ;p)(2^)qg’ {k + p + q) 

■^ang™ = J d^pd^qd^kPi{q-p)Pi'{k-p){‘ 2 TrfS^^\k + p + ^ 

= [ d^pdflqdflkPiiq ‘ p)Pi'ik ‘ q)( 27 T)^ 6 ^^\kp 

(55) 


Note that the second and third integrals above are really 
the same under p ^ q. We term the first integral above the 
symmetric integral and the second and third asymmetric. 
Figure [^explains these equations and their symmetries di- 
agrammatically to illustrate the underlying structure of the 
covariance calculation up to this point. 

To evaluate these angular integrals, we write the Dirac 
delta as the Fourier transform of unity, 

( 27 r)® < 5 ® {k+p + ^ = J d^r ( 56 ) 

expand each exponential in spherical harmonics using 
AWH13 equation (16.61), and perform the angular integral 
over dflr- Defining 

'Piihi3ik,P,q) = J r‘^drji^{kr)ji^{pr)ji^{qr), 

Phhh = (57) 

pj I (2/i + 1)(2/2 + 1)(2/3 + 1) 

= Y-4^-, 

we obtain 


{2TrfS^^\k+p + ^ = 

Plihh^hhh'^hhh (k,p, q) 



h 

l3 

f ll h 

h \ 

1 0 

0 

0 j 

1 mi m2 

m3 J 


^ ^imi (^)^2'^2 (^)^3'^3 (^) ■ 


This is equivalent to Mehrem (2002) equation (5.1) if the 3j- 
symbols above are translated to Clebsch-Gordan symbols. 

Inserting equation (58) into equation (55) and then ex¬ 
panding the Legendre polynomials in equation (55) into 
spherical harmonics using the spherical harmonic addition 
theorem §, we now simply have integrals over products of 
three spherical harmonics, which can be done analytically 
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Figure 4. This shows the symmetries of equation ( |54| , which in 
turn derive from the structure of equation ( |50| ); one can directly 
compare the arguments of the exponentials in this latter with the 
diagram. The leftmost triangle represents the term Ji (q, p)Pi (q-p) 
outside the curly brackets in equation ( |54| ), showing also the ra¬ 
dial arguments implicit in the , and the six triangles inside the 
curly brackets above represent the six terms in the curly brackets. 
The Legendre polynomials are always evaluated about a particu¬ 
lar vertex, as shown in the diagram, and the real-space variables 
match to Fourier-space variables differently in each triangle (see 
equation ( |50| )). One can see from above that each pair of tri¬ 
angles, or pair of terms in curly brackets in equation ( |54| ), has 
switch symmetry ^ f^. These are rotation symmetries about 
the vertex between and in each pair. One can also see that 
if we switch fi ^ r 2 and q ^ p, the leftmost triangle is symmet¬ 
ric. The topmost pair will also be symmetric under this switch 
as well, which is why it gives rise to two symmetric projection 
integrals, but the middle and bottom pairs will not be, which is 
why they give rise to four asymmetric projection integrals (see 
equation ( |62| )). 


with 3j-symbols. The result can then be simplified by ex¬ 
plicitly evaluating some of the 3j-symbols (using NIST Digi¬ 
tal Library of Mathematical Functions (DLMF) 34.3.1) and 
summing over all of the spin angular momenta (using NIST 
DLMF 34.3.10 and 34.3.18). For the symmetric integral we 
find 




(47r)^^(-iy= (2/2 + 1) 

h 


i' 

0 


^^2^20 (P? 0.-) ^) 


(59) 

where we have separated k with a semicolon because it is 
the only argument that does not appear in the Legendre 
polynomials in the integral. A simple case to check is setting 
= 0 and /c = 0 in equation (55). Then q = —p so by direct 
computation 


l™ = I287r®(-1) 

fe-s-o \ ^ 


(p-9) 


(60) 


where we used Pi{—1) = (—1)^- In equation (59), I' = 0 
sets I 2 = I and the 3j-symbol’s square is 1/(2/ + 1). Us¬ 
ing the orthogonality relation for spherical Bessel functions, 
7liio{p,q,0) = 7rS^\p — q)/{2q‘^)] inserting this in equation 


(59) and simplifying yields agreement with the direct com¬ 
putation. 

For the asymmetric integral we hnd 


xJ2i‘2l2 + l)( Q Q qU (-l)(*+*'+*=>/"7^^,^^(p,9,fc), 
h ^ ^ 

(61) 


where now p is separated by a semicolon because it appeared 
in two Legendre polynomials in the integrand. Note that for 
/^ = 0, the symmetric and asymmetric integrals of equation 
(55) are equal, so equation (61) should reduce to equation 
(59) in this limit, as can be verihed by noting /' = 0 implies 
/ = /2. 

Thus 


Cov„, = / I ^^P(p)P(g)P(fc)(2/ + l)(2/' + l)(-l)*+*' 

X Ji iq, P) { (g, (g, P-, k) + J:, (p, {p, g; k) 

+ Jv ik,p)i:ZT(P’ fc) + (P’ k)lZT(P’ fc) 

+ JHk,q)i:Z:wiP,P,k) + JHqT)lZT(P’P’’^')}- ( 62 ) 


We now interchange the order of integration so that 
the integrals over and k are done hrst, since they are 
separable, and the linking integral over r implied hy IZi ^121 s 
is done last. We also make the sum over I 2 explicit and do 
it after evaluating the q^p and k integrals. Finally we dehne 

^Pik)ji{kn)jiikr) (63) 

and 

^P{k)ji{kn)jt,ikr[)ji,{kr). (64) 
In terms of these functions. 


Cov;;/(ri,r2;ri,r2) = -4(2/ + 1)(2/' + !)(-!)*+* 


X 


{(-l)'4o(r) 


fhW {r; ri, r[)fi^w (r; r 2 ,4) 


+ fi2w{i'\r2,ri)fi^u>{r;ri,r2) 




X 

+ 

+ 

+ 


fair-, (r; r 2 , 4) 

fu{r; ri)fvi> (r; rAfi^iv p-, r 2 , ri) 
fu (r; r 2 )fi'i' (r; rAfi^w (r; ri, ri) 

III (r; r 2 )fi'i' (r; rAfhii' (r; ri, ri) 


(65) 


One can see that this is symmetric under switching ri F4 r 2 
and r'l F4 r 2 , as expected. We have thus shown how to re¬ 
duce the covariance of our multipole decomposition to a sum 
of 2-D integrals. This is a signihcant computational beneht: 
the fu and fi^w can be pre-computed once to give all the 
terms in the sum above, and then integrated over dr. Fur¬ 
ther, since the problem is now 2-D one can simply evaluate 
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the integrals using a grid and avoid appealing to more com¬ 
plicated higher-dimensional integration techniques. 

In closing, we note that for P oc 1/k (i.e. oc l/r^), fu 
and fi^ii/ can be computed analytically. We find 

_p ^ \ 1 z -1-2 r(/ +1) 

x^(' + i’y^ + y(^) ) 

using ji(x) = ^/7T/(2x)Ji^i/2{x) and Gradshteyn & Ryzhik 
(2007) equation (6.512.1). F is the hypergeometric function 
and we assume r < ri; the result for r > ri is given by 
switching ri o r. fi^u' can be computed using techniques 
outlined in Fabrikant (2013) and is given by his equation 
(9); since the expression is rather long we do not reproduce 
it here. We mention this since one could imagine scenarios 
in which high speed was desirable for computing the covari¬ 
ance, such that these approximate forms might suffice. Fi¬ 
nally, to incorporate shot noise in the covariance, one takes 
P{k) P{k) + l/n, n the survey number density. This will 
introduce a cross term where one of the fu or fi^iu in each 
pair in equation ( |65| ) no longer involves the power spectrum, 
and also a term in 1/n^ where both functions in each pair 
do not. The required integrals are also analytic: 

^ (67) 

while fi^iu is rather longer and given by Mehrem (2002) 
equation (5.14), assuming a much simpler form (his equation 
(5.15)) if/,/', or h = 0. 


6.4 Radial binning for the covariance 


The above calculation used exact values for ri, r 2 , r'l and 
r 2 , but we can easily integrate over bins in radius. We simply 
integrate the fi^u' and fu functions defined above over bins, 
equivalent to replacing ji{kri) and ji{kr'i) with their bin- 
averaged values. We define 

^Pik)Mk-,n)Mkr) (68) 

and 

/ h‘^rjh 

(fc; ri)ji> (fc; ri)ji^ (fcr) (69) 


with 


ji{k;ri) 


f vfdu ji{ku)^(u] Ti) 
f u^du 4>(i^; Ti) 


(70) 


where u is a dummy variable and we recall that ^(u;ri) is 
the binning function ensuring u is in the bin Vi (see Section 


2.3). 


6.5 Covariance results 

We display the binned reduced covariance. 

Red Covn'(ri,r 2 ;ri,r 2 ) = 

_ CoYu'{ri,r2]r[,r2) _ 

Co\u (r'l, r' 2 ; r [, r^)Covz/z/ (r [, ; r [, r^) 


for fixed ri, r 2 and a number of ri, r 2 and //' combinations in 
Figure]^ One can see clear features when (ri, r 2 ) = (^i, '^ 2 ), 
especially when / = /' as well (e.g. the 11, 22, 33, and 44 
panels). The computation was done in 8 Mpc/h bins but we 
display with an interpolated color scheme because the un¬ 
derlying radial variables are continuous, in contrast to the 
multipoles / and /'. We used the linear-theory matter power 
spectrum from GAME (Lewis 2000) and checked conver¬ 
gence of the integrals by varying the endpoints and spacing 
of the grids in r and k we used[^ For the spherical Bessel 
functions we used high-order Taylor series for small values 
of the arguments, with the change-over point to the series 
depending on the order /, and cross-checked with direct com¬ 
putation using scipy’s built-in functions. For the ji (equation 
( [TQI ) we used analytical results, cross-checked with numeri¬ 
cal integrations of the ji. 

In Figure we show the binned covariance when ri = 
r'l and r 2 = r 2 versus all / and /', and for a number of choices 
of r'l and r' 2 ^ indicated in the upper left of each panel. Notice 
that the strongest covariance is, as one might expect, when 
/ = /' as well, along the diagonal. We display with no color 
interpolation because the multipoles are discrete. We show 
larger radial bins than the LasDamas mock results contain 
because these will be relevant for the Baryon Acoustic Os¬ 
cillation (BAO) scale analysis planned for future work. 


7 MOCK DATA RESULTS 
7.1 Full results 

We present the results of running our algorithm on the pub¬ 
licly available LasDamas mock catalogs for the SDSS-II DR7 
in both real and redshift space We used 15 radial bins with 
Ar = 6 Mpc//i. We show first the results at each multipole 
versus the two triangle side lengths ri and r 2 in Figure 
This shows that the largest amplitude contribution to the 
3PGF, especially for triangles well away from the diagonal, 
is / = 2. This is what we expect from SE15, Figure 9, third 
row, leftmost panel, showing the perturbation theory results 
and focusing on the linear bias hi , which dominates the non¬ 
linear bias 52. In / = 1, there is a hint of a large-scale decre¬ 
ment, to be compared with the slight feature close to the 
diagonal around 130 Mpc(= 90 Mpc/h) in SE15 Eigure 9, 
second row, leftmost panel. 

As in SE15, / = 2 and / = 3 look similar but with / = 2 
having higher amplitude away from the diagonal. Eor / ^ 3, 
the panels all begin to look the same, agreeing with our 
expectation from SE15 Eigure 9. This is because these higher 
multipoles, in particular near the diagonal, are dominated 
by a small population of squeezed triangles where two sides 
are equal (e.g. ri and r 2 ) and the third side nears zero. In the 
hierarchical ansatz for the 3PGE, one has (/ ~ l/(rir 2 )^ + 
l/(7"2T’3)^ + l/(r 3 ri)^, so we expect the amplitude to become 
very large as any side approaches zero. 

Also discussed in SE15 is another reason for the similar¬ 
ity of the / ^ 3 panels: before cyclic summing over vertices of 


^ Here we used the redshift zero power spectrum from a flat 
AC DM cosmology with D^/i^ = 0.0226, Qch‘^ = 0.112, Ug = 
0.96, and ag = 0.821. 

^ http://Iss.phy.Vanderbilt.edu/lasdamas/mocks/ 


© 0000 RAS, MNRAS 000 , 000-000 












An 0{N‘^) Three-Point Function 13 


Reduced 3PCF covariance 
r[ =96—104, r 2 =64—72 Mpc//i 
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Figure 5. The reduced covariance (equation ( |7l| )) for a number of multipole combinations. The r'^ and bins are fixed, and the ri and 
r 2 bins are the horizontal and vertical axes of the plot; the ll' combination is indicated in the upper left of each panel. We have chosen 
an (r[^r 2 ) bin combination for relevance to the BAO scale 100 Mpc/h) while avoiding the squeezed limit (hence away from r'^). 
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Diagonal of 3PCF covariance 



8.0e-10 


4.0e-10 


0.06+00 
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Figure 6. The diagonal of the binned covariance, i.e. equation ( |65| ) binned and with (ri,r 2 ) = The (ri,r 2 ) bins are indicated 

in the up per left of each panel, and the horizontal and vertical axes show the multipoles I and I'. We have used a survey volume V in 
equation ( |65| > of 1 (Gpc/h)^, roughly that of the SDSS DR12, to normalize. This is about 7 times smaller than the total volume of mock 
catalogs used here; thus l/\/7 times the square-root of the diagonals above gives a rough estimate of the error bars we might expect as 
of order 10“^, or about 5% of the 3PCF (comparing with Figure]^. 
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the triangle, the leading order pre-cyclic perturbation the¬ 
ory 3PCF only has structure for / = 0,1, 2. In particular, at 
leading (fourth) order, the 3PCF receives one contribution 
from the second-order density field, 5^‘^\ which is in turn 
calculated by integrating a kernel against the linear 

density field (Goroff et ah 1986; Jain & Bertschinger 1994; 
Bernardeau et ah 2002). This kernel has only / = 0,1, and 2 
terms. If one chooses the second-order density point to be at 
the origin, the 3PCF therefore has multipole structure only 
to / = 2. In reality we do not know which point contributes 
^(2) 

, SO we must cyclically sum around the triangle and can¬ 
not choose at the origin. This cyclic summing generates 
additional angular structure, but it just stems from the ge¬ 
ometric effect of writing a simple I = 0,1, and 2-only multi¬ 
pole expansion with argument e.g. fi • fa = cos O 13 in terms 
of ri ‘ r 2 — cos 0 i 2 - 

We note that the I = 0,1, and 2 terms that enter pre- 
cyclically have a physical meaning. is formed by sum¬ 
ming two mode-coupling kernels a and 13 (Bernardeau et al. 
2002 equations (39) and (156)). These in turn come from 
solving respectively the full continuity equation and the Eu¬ 
ler equation (compare Bernardeau et al. 2002 equations (16) 
and (17) with their equations (37) and (38)). a produces all 
of the I = 0 and 5/7 of the I = 1 terms in F^‘^\ The I = 0 
contribution is from the product of the velocity divergence 
and the density, while the I = 1 contribution is from gradi¬ 
ents of the density field parallel to the velocity. Meanwhile, 
{3 generates the remaining 2/7 of the I = 1 term and all of 
the I = 2 term in these stem from gradients of the 

velocity divergence parallel to the velocity. 

Figures and show that the full 3PCF of the data 
can be reconstructed well from only a few multipoles. Fig¬ 
ure reconstructs the 3PCF from coefficients Q, up to and 
including the I indicated in the legend, for a particular trian¬ 
gle configuration with n = 70 Mpc//i, r 2 = 40 Mpc/h. One 
recovers an accurate shape versus cos ^12 even using only 
multipoles up to / = 6, and that adding in / = 5 — 8 and fi¬ 
nally / = 5 —10 changes the shape very little. In Figurej^ we 
illustrate the same idea for three different triangle configu¬ 
rations: the higher multipoles fall off relative to ^o, meaning 
they contribute less to reconstructing the full 3PCF. This 
plot likely is conservative in that it makes the effect of the 
higher multipoles appear larger than it is; the plot shows 
the ratio of each higher multipole to Co, but the change in 
the 3PCF produced by adding in a higher multipole is ac¬ 
tually roughly the ratio of the multipole to the sum of all 
the lower multipoles. Since, in detail, Legendre polynomial 
weights also enter, one might consider an angle-averaged ver¬ 
sion of this ratio. However since Figure [^effectively already 
shows the unimportance of the highest multipoles, we have 
in Figure [^ just chosen to show |0/Co| because it offers more 
granular information. 


7.2 Compressing the data 

SE15 presented a compression scheme for the multipole 
moments of the 3PCE. This was designed to avoid the 
squeezed limit where two galaxies are so nearby that per¬ 
turbation theory is invalid and also to reduce the dimension 
of the covariance matrix required for parameter fitting. This 
approach integrated each multipole moment over r 2 from 
ri/3 < r 2 < 2ri/3 at each value of ri. 


Summed multipoles; 



Figure 8. This shows reconstruction of the full 3PCF (for the 
LasDamas real space mocks) from its multipoles using coefficients 
up to and including the I indicated in the legend, as described in 
Section The reconstruction converges even for modest 1. The 
I = 8 points lie essentially directly under those for I = 10. 


In the current work the data is binned coarsely enough 
in both n and r 2 that this approach must be adapted 
slightly. We simply choose to, for a given bin ri, sum over 
all bins with r 2 G S{ri). S{ri) is the set of all bins in r 2 
where r 2 is greater than 3Ar and less than n — 3Ar|^This 
assures that the minimum value of r 2 is 18 Mpc//i and that 
the minimum difference between n and r 2 is also 18 Mpc//i, 
meaning by the Triangle Inequality that ra ^ 18 Mpe/h. 
This avoids the squeezed limit while reducing the dimension 
of the problem. Mathematically, the compression is defined 
here as 


^ Er2eS(ri)C!(n,r2)Ar(r2) 


(72) 


where bar denotes “binned”, superscript “c” denotes “com¬ 
pression”, and Ci(ri,r 2 ) is the l^^ binned 3PCE multipole 
(see Section [ 2 ^. Ay(r 2 ) is the volume of bin r 2 . The de¬ 
nominator is for normalization. 

In Eigures [To] and m we show the results of this com¬ 
pression. We also compressed the leading (fourth) order per¬ 
turbation theory predictions, calculated as outlined in SE15, 
and show them for comparison. The theory requires linear 
(5i) and non-linear (^ 2 ) bias parameters as an input; we use 
a least-squares fit with points weighted by the inverse com¬ 
pressed variance. This latter is computed from the scatter 
between mocks and ignores noise in the random catalog used 
for edge correction, which due to the large number of ran¬ 
doms is negligible. Our mocks constitute a volume of order 


Note that if one used different bin sizes Ar one might wish to 
select a different multiple of Ar in defining S. 
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LasDamas Realspace 3PCF^ ^ 

^^-^ 1 —^ ^ 6.0e-04 
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Figure 7. Binned multipole coefficients Q of the 3PCF as defined in equations Q and ( |l2|). These are the result of applying our 
algorithm to the LasDamas real space mock catalogs for the SDSS DR7, as described in Section]^ and with edge correction as described 
in Section^ The horizontal and vertical axes are ri and r 2 in Mpc/h. As in SE15, we have weighted by the volume in each spherical shell 
r^r|/(100 Mpc/h)^; this is to amplify the finer features. Note that on the diagonal, we expect the 3PCF to be dominated by squeezed 
triangles for which perturbation theory is not valid, so we have not computed the diagonal. For this reason we need not include the 
additional exponential suppression of the diagonal used in SE15. 0 0000 RAS, MNRAS 000, 000-000 
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Figure 9. This shows the ratios of / > 0 coefficients to (^o for sev¬ 
eral triangle configurations (again using the LasDamas real space 
mocks). The decline of the higher multipoles with I indicates that 
not many multipoles are needed for accurately reconstructing the 
full 3PCF. This is especially true for the largest scale triangle we 
show, which is also the least likely to be altered by non-linear ef¬ 
fects. The relative magnitudes of the higher multipoles here may 
seem large when recalling from Figure that the reconstruction 
appears well converged by / = 6; but note that a given multipole’s 
contribution to the reconstruction is roughly its ratio to the sum 
of all the lower multipoles, not just to Co; this reduces the im¬ 
portance of the higher multipoles. Finally, the strength of / = 2 
shows the quadratic or U-shaped behavior of the 3PCF tradition¬ 
ally associated with gravitational growth of structure (see also 
Figure]^. Gravity generates gradients of the density and veloc¬ 
ity divergence mostly parallel to the velocity, in turn enhancing 
roughly collinear structures with 6 near 0 or tt (Bernardeau et al. 
2002 ). 


7 times that used for the theoretical covariance calculation 
here, so, as explained in Figure we might expect error 
bars on the compressions of order 5%. This is indeed what 
we find. We offer the caveat that a full, rigorously correct fit 
of theory to observation would require inversion of the full 
covariance matrix. We leave this for future work; here our 
goal is simply to indicate that the results of our algorithm 
roughly agree with perturbation theory predictions. 

Using the simple procedure above, the results are well- 
fit with hi — 1.90 and ^2 = 0.93; note that to compute the 
theory predictions we matched the LasDamas cosmology]^ 
There is some deviation noticeable at large scales in / = 0 
(about 3cr) and I = 1 (about 2cr) , with nearly all the other 
multipoles deviating only within the error bars or at most in 
a few cases just slightly outside them. The larger deviations 
in / = 0 and I = 1 are likely because non-linear corrections 



q [Mpc/ft.] I 


Figure 12. The left panel shows the ratio of redshift space to real 
space results for the LasDamas mocks at each multipole and in 
each radial bin whose compression is non-zero (ri ^ 45 Mpc/h). 
The right panel shows the radial average at each multipole. While 
there is some radial scale dependence in the left panel, it is modest 
for all but / = 8 — 10. Thus averaging over the radial dependence 
does not lose much information for the lower multipoles. The 
averages (right panel) are similar for all but I = 8 — 10, and even 
these differ by less than a factor of 2 from the averages of the 
lower multipoles. 


to the perturbation theory results cannot be neglected. In 
particular, in / = 0 we expect non-linear evolution might 
smooth structure on smaller scales, making the slope of the 
perturbation theory compression shallower and allowing a 
better global fit to the I = 0 mock results. 

Importantly, the error bars become much larger for I ^ 
5 as compared to those for / < 5. This suggests when doing a 
full parameter fit using the covariance matrix, one might not 
gain much by including these higher multipoles. One might 
choose simply to drop these modes to reduce the dimension 
of the covariance matrix to be inverted. 

Figures and m also show that the redshift space 
results at each multipole appear to be roughly a constant 
rescaling of the real space results, with a constant that only 
weakly depends on the multipole. To illustrate this we show 
the ratio of redshift space to real space results in each ra¬ 
dial bin at each multipole (Figure left panel) and the 
radially-averaged ratio versus multipole (Figure |12[ right 
panel). More detailed discussion is in the caption to this 
Figure; the key point is that for / < 8, there is little ra¬ 
dial dependence to the rescaling factor and also little mul- 
tipole dependence. Both dependences are more pronounced 
for / = 8 — 10; we suspect this is because these higher mul- 
tipoles are dominated, even in the compression, by a small 
subset of relatively squeezed triangles that are more strongly 
affected by RSD. This issue might merit further attention 
in subsequent work. 


^ LasDamas cosmology given at: http://lss.phy.vanderbilt. 
edu/lasdamas/simulations .html ag = 0.8 there is quoted at 
z = 0; Dm = 0.25, Da = 0.75, rig = 1. The LasDamas mocks are 
at z = 0.3, so when normalizing the power spectrum we should 
use (Tsiz = 0)[D(0.3)/D(0)], with D the linear growth factor (e.g. 
Mo van den Bosch & White (2010) equations (4.75), (4.76), and 
(3.77); Carroll et al. (1992)). 


7.3 Compressing the covariance 

Given that we now have compressed data, we must also ap¬ 
ply our compression scheme as in the previous section to 
the binned covariance of Section 12.31 We denote the com¬ 
pressed, binned covariance as Cov^^/(ri; rj), noting that the 
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1=0-4 compressions, =1.90, 63 =0.93 



Ti [Mpc/h] 


Figure 10. The results of compressing the binned multipole moments as described in Section [7.2| The perturbation theory points are 
predictions for the 3PCF at lowest (fourth) order, with linear bias bi and non-linear bias bi given in the title (see e.g. SE15 equations 
(1), (2), and Sections 6 and 7; note our non-linear bias does not take a factor of 1/2 when it multiplies the matter density field’s square). 
Intrinsically the 3PCF is of order with ^ the 2-point correlation functon, so on large scales we expect it to be of order 10“"^. We 
therefore multiply the axis labels by this factor for compactness. Importantly, notice that the redshift space mock results are roughly 
just a rescaling of the real space mock results by a radius-independent constant that only weakly ( 
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1=5-10 compressions, =1.90, 63 =0.93 



Vi [Mpc/h] Ti [Mpc/h] 


Figure 11. Same as Figure [To] but for the higher multipoles, 1 = 5 — 10. Note the lower amplitude of these as compared with especially 
the lowest multipoles in Figure [T^ and the larger errorbars. 
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two superscript “c”s denote that we compress over r 2 and r 2 , 
leaving the quantity a function only of ri and r'l. 

It would be computationally intensive to compute the 
binned covariance using equation (65) with the fu and fi^u^ 


replaced by equations (68)-(70) and then compress, and a 


faster approach is available. This is to compress fu and 
fi^iu as necessary first and from them obtain the compressed 
binned covariance, fu need only be compressed at most once 
(if its argument is r 2 or r^), but fi^u' may be compressed 
once or twice depending on if one or both of its arguments 
have subscript 2. We thus define three functions, where “cc” 
again denotes a double compression: 


mr-,n) = 




ft 2 llir;ri,r[) = 


Er2SS(ri)^^(^2) 

Er^sS(ri) fhii' {r-, ri , ra) AV(rj) 


Er^SS(r'j) ^^(■'’2) 

Er2gS(ri)Er'gS(r;)4tt'(^;^2,r^)Ay(r^)Ay(r2) 

Er2gS(ri) 


(73) 


Note that in the second line above, the fi^u' being com¬ 
pressed is a function ri and r 2 , and hence need only be 
compressed once—it is not compressed over ri. However in 
the third line above, the fi^u' being compressed depends on 
r 2 and r 2 and so must be compressed twice. Making these 


replacements as appropriate in equation (65) yields the com¬ 


pressed, binned covariance. This shows that the framework 
of compression can be easily generalized to the covariance. 


8 CONCLUSIONS 

We have presented a novel algorithm to compute the multi¬ 
pole moments of the 3PCF. It is especially apt for large cos¬ 
mological datasets such as SDSS and upcoming surveys like 
Euclid, Large Synoptic Survey Telescope (LSST), and Dark 
Energy Spectroscopic Instrument (DESI), which will have 
tens of millions to billions of objects (Jain et al. 2015). Eor 
these datasets, an approach that scales with would be 
wholly infeasible. We have shown that our algorithm scales 
as , handles edge correction easily, and permits computa¬ 
tion of the 3PCE of a large dataset quickly even with modest 
computing resources. We have also computed the covariance 
matrix of this decomposition in the Gaussian random field 
limit. Einally, we have developed the compression scheme 
first presented in SE15 and shown its application both to 
the data and to the covariance matrix. This compression 
scheme offers a compelling way to visualize the results of 
the algorithm that loses little information, in contrast to the 
plots of the 3PCE or reduced 3PCE versus opening angle 0 
for particular triangle configurations that previous literature 
supplies. 

The algorithm presented here is unique in that it fun¬ 
damentally reduces the scaling of the 3PCE measurement to 
that of the two-point function, while remaining exact in an¬ 
gle. This did not have to be the case. Eormally, for a complete 
representation of the 3PCE, one needs an infinite number of 
multipoles L However, because the physics generating the 
3PCE does not have a great deal of angular structure, in 


practice a finite, modest number of multipoles suffices. Eur- 
thermore, we have shown that in our LasDamas test case, 
the 3PCE is already well-reconstructed by / = 6 (Eigurej^. 
Since our algorithm fundamentally requires pair-counting, 
using a fast Eourier transform (EET) for this step may in 
some cases offer an additional acceleration; we present this 
in Slepian & Eisenstein 2015c. 

Einally, one might worry that jagged survey boundaries 
could easily introduce high multipoles into the measured 
3PCE. But we have shown that the coefficients required for 
the edge correction, at least for our LasDamas test case, 
fall off quickly enough that one need only measure a few 
multipoles of the randoms for an accurate solution (Eigures 
and . 

The 3PCE contains important information on the non- 
Gaussianity of large-scale structure (LSS) due to growth 
under gravity and also perhaps that remaining from pri¬ 
mordial non-Gaussianity. With measurements of only the 
2 -point function, the amplitude of clustering (e.g. erg) and 
the linear bias are degenerate. However, the 3PGE is sensi¬ 
tive to a different power of the linear bias than the 2 -point 
function (cube versus square), and so measuring it exposes 
a raw factor of the bias and helps break this degeneracy. 
As for primordial non-Gaussianity, while the GMB has been 
the dominant constraint up to now (see Ade et al. (2015)), 
it is expected that even maximally improved GMB measure¬ 
ments can only enhance the GMB constraint by a factor of 
a few. Thus LSS will become a vital complementary probe. 
Generically, inflation must couple to ordinary matter so as 
to produce it during reheating, and this coupling produces 
some level of non-Gaussianity (Desjacques & Seljak 2010, 
for a recent review). Thus the 3PGE can be used to probe 
the dynamics of inflation in principle—and perhaps, soon, 
in practice. 

The 3PGE also contains information on redshift space 
distortions. It should be emphasized that in the current 
work, the multipole moments are averaged over rotations 
of the triangle configurations, and so we lose any informa¬ 
tion about anisotropy. However, our algorithm can easily be 
adapted to retain the full, unaveraged information (aims) 
around each possible origin. This would allow tracing the 
anisotropy RSD induce. Preliminary calculations indicate 
that RSD introduce couplings between multipoles I 7 ^ I' 
which are absent without RSD. These couplings have se¬ 
lection rules due to the underlying symmetries under ro¬ 
tation about the line of sight and parity flips. There will 
also be m dependence induced by the preferred direction 
defined by the line-of-sight. We therefore expect that the off- 
diagonal terms in a tensor of spherical harmonic coefficients 
will have structure that can be used to probe RSD-induced 
anisotropies. It is already important that the spherical har¬ 
monics, with the introduction oil' ^ I and m, offer a natural 
5-D basis for redshift-space measurements. Work on these 
questions from both analytic and numerical perspectives is 
underway. 

We plan to apply our algorithm and analysis approach 
to SDSS DR 12 , with the goals of assessing the presence of 
BAG features, measuring the linear and non-linear bias, and 
constraining the baryon-dark matter relative velocity (Tseli- 
akhovich & Hirata 2010; Yoo, Dalai & Seljak 2011; Yoo & 
Seljak 2013; SE15). Previous literature has found a BAG 
feature in the reduced 3PGE with ^ the 2PGE (Gaz- 
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tanaga et al. 2009). That work used only one triangle config¬ 
uration, ri = 33±5.5 Mpc//i and r 2 = 88±5.5 Mpc//i. With 
the additional signal-to-noise the large number of galaxies 
in SDSS DR12 offers, as well as our algorithm’s ability to 
consider all triangle configurations quickly, this problem is 
ripe for revisiting. Furthermore, SE15 suggests that the mul- 
tipole decomposition clearly isolates a strong BAO feature, 
especially in the I = 1 multipole. This is also a particu¬ 
larly informative multipole for the relative velocity, further 
discussed in SE15. That work additionally shows that, in 
principle, the multipole decomposition can clearly separate 
the effects of linear and non-linear bias—significant because, 
as noted above, the 3PCE has traditionally been an impor¬ 
tant tool for constraining these parameters. Einally, the sig¬ 
nificant speed advantage of our algorithm will permit much 
finer and much faster calibration of any 3PCE measurements 
against large cosmological simulations. Such improved cali¬ 
bration should greatly enhance the leverage of the 3PCE as 
a fundamental probe of large-scale structure. 
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